function [ out ] = fun_d2sigma_new(theta,argument,m,combos_all,ind1,ind2)

v = combos_all - 1;

v_minus1_minus1 = v;
v_minus1_minus1(:,[ind1 ind2]) = v_minus1_minus1(:,[ind1 ind2])-1;

v_plus1_plus1 = v;
v_plus1_plus1(:,[ind1 ind2]) = v_plus1_plus1(:,[ind1 ind2])+1;

v_minus1_0 = v;
v_minus1_0(:,ind1) = v_minus1_0(:,ind1)-1;

v_plus1_0 = v;
v_plus1_0(:,ind1) = v_plus1_0(:,ind1)+1;

v_0_minus1 = v;
v_0_minus1(:,ind2) = v_0_minus1(:,ind2)-1;

v_0_plus1 = v;
v_0_plus1(:,ind2) = v_0_plus1(:,ind2)+1;

v_minus1_plus1 = v;
v_minus1_plus1(:,ind1) = v_minus1_plus1(:,ind1)-1;
v_minus1_plus1(:,ind2) = v_minus1_plus1(:,ind2)+1;

v_plus1_minus1 = v;
v_plus1_minus1(:,ind1) = v_plus1_minus1(:,ind1)+1;
v_plus1_minus1(:,ind2) = v_plus1_minus1(:,ind2)-1;

v_0_0 = v;

B_der = prod(BernsteinElement_2(v_minus1_minus1,m,argument'),3).*(m(ind1)-(v(:,ind1)-1)).*(m(ind2)-(v(:,ind2)-1)) +...
         prod(BernsteinElement_2(v_0_minus1,m,argument'),3).*(2*v(:,ind1)-m(ind1)).*(m(ind2)-(v(:,ind2)-1)) +...
         -prod(BernsteinElement_2(v_plus1_minus1,m,argument'),3).*(v(:,ind1)+1).*(m(ind2)-(v(:,ind2)-1)) +...      
        prod(BernsteinElement_2(v_minus1_0,m,argument'),3).*(m(ind1)-(v(:,ind1)-1)).*(2*v(:,ind2)-m(ind2)) +...
         prod(BernsteinElement_2(v_0_0,m,argument'),3).*(2*v(:,ind1)-m(ind1)).*(2*v(:,ind2)-m(ind2)) +...
         -prod(BernsteinElement_2(v_plus1_0,m,argument'),3).*(v(:,ind1)+1).*(2*v(:,ind2)-m(ind2)) +...
         -prod(BernsteinElement_2(v_minus1_plus1,m,argument'),3).*(m(ind1)-(v(:,ind1)-1)).*(v(:,ind2)+1) +...
          -prod(BernsteinElement_2(v_0_plus1,m,argument'),3).*(2*v(:,ind1)-m(ind1)).*(v(:,ind2)+1) +...
         +prod(BernsteinElement_2(v_plus1_plus1,m,argument'),3).*(v(:,ind1)+1).*(v(:,ind2)+1);

out = B_der'*theta;

end